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In this paper we study the electrostatic properties of 'Janus' spheres with unequal charge densities on both 
hemispheres. We introduce a method to compare primitive-model Monte Carlo simulations of the ionic double 
layer with predictions of (mean-field) nonlinear Poisson-Boltzmann theory. We also derive practical DLVO- 
like expressions that describe the Janus-particle pair interactions by mean-field theory. Using a large set of 
parameters, we are able to probe the range of validity of the Poisson-Boltzmann approximation, and thus 
of DLVO-likc theories, for such particles. For homogeneously charged spheres this range corresponds well to 
the range that was predicted by field-theoretical studies of homogeneously charged flat surfaces. Moreover, 
we find similar ranges for colloids with a Janus-type charge distribution. The techniques and parameters we 
introduce show promise for future studies of an even wider class of charged-patterned particles. 
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I. INTRODUCTION 

Electrostatic interactions in suspensions of charged col- 
loids are of paramount importance to the structure and 
phase behaviour of such systems 1 - - — . The Dcrjaguin Lan- 
dau Verwey Overbeek (DLVO) 13 ' 14 theory is the most 
well-known theory by which the pair interactions between 
screened charged particles can be described. However, 
DLVO theory only describes the interactions between 
homogeneously charged spherical particles and is there- 
fore a monopole theory. The rapidly-growing zoo 15 ' 16 of 
new colloidal particles demands similar theories that are 
equipped to deal with anisotropy in shape and size. 

For classical, spherically charged, particles, DLVO the- 
ory has proven to be a powerful means to describe 
the electrostatic properties of systems; mostly for high- 
polarity solvents, low surface charge, and high ionic 
strength. The nonlinear Poisson-Boltzmann (PB) ap- 
proach^— extends this range, although analytic re- 
sults are only possible for a limited number of systems. 
PB theory is based on a mean-field approximation that 
ignores ion-ion correlations, which is justified only for 
sufficiently high temperatures, not too apolar media, 
and monovalent ions at reasonable concentrations. In 
regimes where ion-ion correlations are important Strong- 
Coupling (SC) theory 2 ^"— may be applied. Several mod- 
ifications to PB theory exist&2^~— , including modifica- 
tions of the traditional PB theory that account for finite- 
size ions. 

In the computational field, both Monte Carlo (MC) 
and Molecular Dynamics (MD) techniques have been de- 
veloped to analyse charged particles suspended in an elec- 
trolyte. In primitive-model simulations, ions are taken 



a 'Elcctronic mail: j.dcgraafl@uu.nl 
^Electronic mail: r.vanroij@uu.nl 



into account explicitly, whilst the solvent is modelled as 
a dielectric continuum 3 - 1 - - — . Ewald Sums 3 -^ usually pro- 
vide the basis for calculating the long-ranged Coulomb 
contribution to the total energy of a system with peri- 
odic boundary conditions. However, employing Ewald 
Sums is computationally expensive and the systems that 
can be studied in the primitive model are therefore typ- 
ically small. When studying charged colloids suspended 
in an electrolyte it is desirable to coarse-grain the sys- 
tem and use the much shorter ranged effective interac- 
tions between the particles instead of accounting for the 
ions explicitly. Most simulation studies therefore con- 
sider systems where the interactions between the col- 
loids can be modelled by a DLVO pair potentia l 12 ' 38 " —. 
We note that faster Ewald-Sums implementations exist, 
e.g., mesh-based approache s 42 ' 43 , but these cannot out- 
perform a course-grained method that does not take the 
ions into account at all in a system of charged colloidal 
particles. 

In this paper we investigate the range of applicabil- 
ity for nonlinear PB theory to accurately describe the 
behaviour of the ion density around charged heteroge- 
neous particles. This allows us to quantify the parameter 
regime for which a (multipole-expandcd) DLVO approx- 
imation may be applied to describe pair interactions in 
coarse-grained simulations, since (the electrostatic part 
of) such DLVO theories can be derived using PB approx- 
imations. Charge-patterned particles have already been 
studied using the DLVO approximation by partitioning 
the surface charge over a finite number of point- Yukawa 
charges with different prcfactors 4 ^ to obtain effective pair 
interactions. For charge-patterned particles this approx- 
imation still generally results in an expensive calculation 
of the pair interaction as a function of the separation and 
orientation. Moreover, the point- Yukawa description in- 
adequately accounts for the hard core of the particle, 
i.e., it implicitly assumes that ions can penetrate the col- 
loid^. In this paper a correct, simple, and therefore com- 
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putationally far more efficient DLVO-multipolc approx- 
imation is derived for the effective interaction between 
two charge-patterned particles. The multipole-based ef- 
fective interactions have an enormous potential for use 
in simulation studies to explore the phase behaviour of 
previously inaccessible systems. 

For this study we focus on particles with a Janus-type 
charge pattern. The term Janus refers to the two-faced 
Roman god of doors and was introduced to describe col- 
loid properties in 1988^. A Janus particle^— consists 
of two opposing parts (usually hemispheres) with differ- 
ent properties for the wetting, charge, chemical function- 
ality, etc. The past decade has seen a marked increase 
in the ability to synthesize such Janus colloids 4 *^— and 
their use in self-assembly experiments. Many interest- 
ing structures have been foun d 44 ' 56 and questions have 
been raised on how to approach simulations of such sys- 
tems. With our study we aim to address some of these 
questions for charged Janus particles in an electrolyte, in 
much the same way as the pioneering simulation studies 
that probed the applicability of the common DLVO/PB 
approximation for homogeneously charged particle a 31 ' 33 . 

In Section [TT] we introduce the methods by which we 
compute the ion density around charged Janus particles: 
primitive-model MC simulations pi A|) and nonlinear PB 
theory (jll B[) . We discuss the results of our investigation 
in Section IIII1 which is divided into four parts. In Sec- 
tion IIII Al we introduce the method, based on Fourier- 
Legendre (FL) mode decomposition, by which we com- 
pare the MC and PB results. This method is applied for 
a homogeneously charged particle in Section [ill Bl where 
we also investigate the relation to the field-theoretical re- 
sults of Refs. [H3 and HH for homogeneously charged flat 
surfaces. In Section llll Cl we extend our results to a Janus 
dipole and show that there is remarkable correspondence 
with the results for a homogeneously charged sphere. We 
consider a particle with a single charged hemisphere in 
Section ITlID I Throughout these sections we give explicit 
recipes for calculating the pair interactions between such 
particles within the PB approximation that we are test- 
ing. The interested reader is referred to the Appendix 
for a derivation of these pair interactions. We discuss 
our findings, comment on the potential synergy between 
simulation methods and theoretical results, and present 
an outlook in Section IIVI 



II. SIMULATIONS AND THEORY 

In the following we consider a system of spherical 
charge-patterned colloids with radius a suspended in an 
electrolyte. The colloid volume fraction is denoted by rj. 
We studied three types of charge distribution for the col- 
loids, (i) A homogeneous surface charge of Ze, with Z > 
the number of charges and e the elementary charge, (ii) 
A perfectly antisymmetric surface charge, with charges 
Ze/2 and —Ze/2 homogeneously distributed over the 
particle's upper and lower hemisphere, respectively, (iii) 



A homogeneously charged upper hemisphere with charge 
Ze and a completely uncharged lower hemisphere. Unless 
stated otherwise Z = 100 throughout this paper. We as- 
sume that there is a perfect dielectric match between the 
colloid, the ions, and the medium to avoid any dielectric 
boundary effects. 



A. Ewald Sums and Monte Carlo Simulations 

To study the systems described above by MC simu- 
lations we turn to the primitive model, for which the 
ions are represented by charged spheres and the solvent 
is treated as a dielectric continuum. To simplify the cal- 
culations we study only one of these particles, which we 
locate at the centre of a volume V cc \\ = 47ra 3 /(3?/). We 
apply periodic boundary conditions to this volume to ac- 
count for the fact that we are in principle interested in a 
system which contains many colloids. The particle's (het- 
erogeneous) surface charge is specified by f 00 charge sites 
distributed over this surface, which can be positively or 
negatively charged, or which do not have charge. These 
charge sites on the colloid are chosen according to the 
optimal packing of 100 points on a sphere 5 ^, to ensure 
that they are spaced as homogeneously as possible. 

The number of free monovalent ions N in the volume 
V^eli is fixed, i.e., we are interested in an average ion 
concentration N/V ce n for the system that we approximate 
by our one-colloid calculation. We only consider systems 
for which a monovalent salt has been added to enhance 
the screening effected by the counter ions to the particle's 
charge. The balance between the number of positive N+ 
and negative N_ ions (N = N + + N- ) is such that the 
volume, and thereby the entire system, is charge neutral. 
For the monopole and charged hemisphere we require Z+ 
N + - N- = and for the Janus dipole N + = 7V_ . 

To sample phase space Monte Carlo (MC) simula- 
tions are performed in the isothcrmal-isochoric (canon- 
ical, NVT) ensemble. We consider a cubic simulation 
box of length L = 50d (77 = (47r/3)(a/L) 3 , V cc n = L 3 ), 
for which we employ periodic boundary conditions. Here 
d is the ion diameter and we assume all ions to be the 
same size. The spherical colloidal particle is located at 
the centre of the box (the origin). The particle's rota- 
tional symmetry axis is chosen parallel to one of the box' 
ribs for the Janus-type charge distributions. 

The ion-ion pair potential is a combination of a 
Coulomb and a hard-core interaction part: 

7./ TT fr r 1 - qiQje2 1 if 00 ' \ r i~ r j\< d > m 
U ^^~ Wlr.-r.lM 0, \t t -T j \>d, (1) 

with Vi and Yj the position of ions i and j with respect to 
the colloid's centre, respectively. The function | • | gives 
the Euclidean norm of a vector, qi = ±1 the sign of the 
i-th ion's charge, eo the permittivity of vacuum, and e 
the relative dielectric constant of the medium, ions, and 
particle. The interaction between a charge (^e on the 
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particle located at r$, with |r,| = a — d/2, and an ion 
with charge qje located at rj is given by 



Wsi(ri,rj-) 



1 



47re e r* 



oo, |r,-| < a + d/2; 
0, |r,-|>a + d/2. 



(2) 

The coupling between periodicity and the (long-range) 
Coulomb interactions, Eqs. ((TJ and (|2]l. is taken into ac- 
count using Ewald Sums with conductive boundary con- 
ditions^^. The total electrostatic energy Uc of a par- 
ticular configuration may be written as 
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(3) 



where summation is over both the ions and the charge 
sites, i.e., N is the total number of charges in the sys- 
tem, both free and fixed; k = (2n/L)l is a Fourier space 
vector, with 1 £ Z 3 ; 7 is the Ewald convergence param- 
eter—; and erfc(-) is the complementary error function. 
One can safely ignore the site-site interactions in Eq. ^ , 
because this gives a constant contribution to the electro- 
static energy Uc- The self-energy term also drops out of 
the energy difference, on which the acceptance criterion 
for the MC trial moves is based^. 

For our simulations we employ the following parame- 
ters, (i) Each run consists of 100,000 MC equilibration 
cycles, where 1 MC cycle is understood to be one trial 
(translation) move per free ion. (ii) This equilibration is 
followed by a production run of 250,000 MC cycles to de- 
termine the ensemble-averaged ion density profiles p± (r) , 
with r the position with respect to the centre of the col- 
loid, (iii) The step size for the ion translational moves 
is in the range [0, 5d] and it is adjusted to yield an ac- 
ceptance ratio of 0.25. (iv) For the Ewald Sums the real 
space cut-off radius for the third term in Eq. ([3]) is set 
to L/2.5, and we use 7 = 0.03 and |1| < 6. This choice 
of Ewald parameters gives a reasonable approximation 
to the value of the electrostatic interaction energy Dou- 
bling and halving the number of cycles for several runs 
showed that the MC parameters give sufficiently equi- 
librated results for most systems. A possible exception 
to the perceived equilibration is deep inside the strong- 
coupling regime, where ion-ion correlations play an im- 
portant role, as we will explain in Section riff Al 



B. The Poisson-Boltzmann Approach 

The spherical particle of radius a in a cubic box 
(L x L x L) models a system with colloid volume fraction 



i] = (47r/3)(a/i) 3 . The equivalent system in PB the- 
ory is described using a spherical Wigner-Seitz (WS) cell 
mode l 19 i 61 ~ — , where the radius of the WS cell is given by 



R = 



3£^ 

47T 



-1/3 



(4) 



and the colloid is located at the centre of the cell. The 
choice of R ensures that the volumes, and therefore the 
average density of colloids/ions is the same as in the cubic 
box of the MC simulations. PB theory is a pplied, in 
accordance with the procedure outlined in Refs. l64T[6a . to 
determine the dimensionless electrostatic potential <j>(r) 
and the associated ion density profiles p±(r) around the 
colloid. 

In our MC simulations the hard-core interaction be- 
tween the ions and the colloid prevent the ions from 
approaching the colloid's centre closer than a distance 
of a + d/2. We therefore assume the same spherical 
hard-core exclusion volume for the point ions in PB the- 
ory. The colloid's surface charge density is given by g(r), 
which is only nonzero when |r| = a; the spatial integral 
over g(r) gives the total colloid charge. The PB equation 
for this system may now be written as 



V 2 <j>(v) = 47rA B g(r) + 







|r| < a + d/2 



k 2 sinh(0(r)) |r| > a + d/2 



(5) 

where k 2 = 87tAb/0 s (such that k 1 is the Debye screening 
length), with p s the (yet unknown) bulk ion density and 
Ab the Bjerrum length 



A B 



(6) 



with fc B Boltzmann's constant and T the temperature. 
We impose the following boundary condition 



V0(r) 



0, 



(7) 



with r = r/|r|, on the edge of the spherical cell to ensure 
that the normal component of the electric field vanishes 
at the boundary, i.e., the WS cell is charge neutral. 

To solve Eq. ([5|) with the above boundary conditions 
the charge density q(r) and the electrostatic potential 
4>(r) are expanded into a complete set of Lcgcndrc poly- 
nomials as 



g(r) = ^2aiS(r - a)P t (x); 

e 

0(r)=X>to^(s). 



(8) 
(9) 



with a 1 and 4>e(r) the surface-charge and potential 
modes, respectively. Here r = |r| and x = r • z, with z 
the orientation of the colloid's rotational symmetry axis. 
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The Pi{-) are £-th order Legendre polynomials, i.e., 



III. IONIC SCREENING OF JANUS PARTICLES 



Po(x) = 1; 
Pi (a;) = x; 

P 2 (x) = X - (3x 2 - 1) ; 
P 3 (x) = i(5x 3 -3x) 



P<(z) 



fc=0 



/ l+fc- 
2 



(10) 

(ii) 

(12) 
(13) 

(14) 



where the expressions between brackets in Eq. Q14p de- 
note binomial coefficients. The nonlinear PB equation 
[Eq. (|SJ)] is likewise expanded using Fourier-Legendre 
(FL) mode decomposition and Taylor expansion of the 
sinh(-) term around the monopole potential 4>o(r). The 
higher-order expansion coefficients contain products of 
Legendre polynomials which effects couplings between 
the various modes. These products must be rewritten as 
a sum of single Legendre polynomials^ to solve for the 
separate modes using an iterative procedure. The mode 
coupling this induces necessitates the analysis of a signif- 
icant number of multipoles even if, for example, only the 
dipole mode is of interest. References [45| and [68| can be 
consulted for more information on the procedure of mode 
expansion to solve the PB equation for hetcrogeneously 
charged colloids. 

It is important to note that the PB theory treats 
the screening ions in the grand-canonical (p,VT) ensem- 
ble. The MC simulations were however performed in the 
canonical ensemble, where the number of ions is fixed, 
to allow for faster exploration of phase space. We fit the 
bulk ion concentration p s in PB theory to ensure that 
the number of positive and negative ions in the WS cell 
corresponds to the number of ions in the MC simulation 
box. We consider this condition, coupled with the fact 
that we study the same colloid volume fraction r\ in both 
approaches, sufficient to justify comparison of the results 
in the two ensembles. The bulk ion concentration is fitted 
according to the criterion 



N± = iV±,p B = J drp ±iPB 



(r), 



(15) 



where the integration is over the region |r| E [a + d/2, R]. 
One of the two equations is redundant, since solving for 
N + is equivalent to solving for N-. The appropriate 
bulk ion concentration p Sl which comes into the right- 
hand side of Eq. (fTS"]) via the dependence of p±,pb(i") = 
p s exp(=F</>(r)) on this concentration, is established using 
an iterative procedure. All PB results presented in this 
paper were obtained on an equidistant radial grid of 2,000 
points for |r| £ [a+d/2, R] by 5-th order Taylor expansion 
of sinh(-) using 6 multipolc modes. 



In this section we describe our results for the compar- 
ison of ion density profiles obtained by MC simulations 
and by PB theory. A total of 99 systems are consid- 
ered for each of the three charge-patterned colloids. We 
use three particle radii a = 5c?, 10c?, and 15c?. For every 
particle radius a, three salt concentrations are studied: 
125, 250, and 375 monovalent cations and anions, re- 
spectively, are added to the counter ions already present 
in the system. This gives N± = 175, 300, and 425 for the 
Janus dipole. For the homogeneously and hemisphcri- 
cally charged particle N + = 125, 250, and 375, when 
N- = 225, 350, and 475, respectively. We consider 11 
Bjerrum lengths A B /c? = 0.01, 0.05, 0.1, 0.25, 0.5, 1.0, 
2.0, 4.0, 6.0, 8.0, and 10.0 for each a and N± combina- 
tion. 

We subdivided the results sections on the basis of the 
surface-charge pattern studied. For each type of charge 
pattern, we consider the range where PB theory accu- 
rately describes the system. We also give a multipole- 
cxpanded DLVO description for the effective pair poten- 
tial between two such objects (for the Janus dipole and 
charged hemisphere) , which may be used in this range to 
model the pair interactions in coarse-grained simulations. 



A. Method of Comparison 

Figure [1] shows an example of our results for a typical 
set of parameters: a = 10c?, As = d, and 250 added an- 
ions and cations, respectively. The azimuthal average of 
the net ionic charge [p + (r) — p-(r)] is shown for the Janus 
dipole (Fig. QJi) to give insight into the shape of the den- 
sity profiles. The multipolc-cxpanded cation-density pro- 
files (Fig. [Tb; p+,i{r) for < t < 3) further illustrate the 
level of correspondence between the MC and PB result 
for this system. Due to the antisymmetry of the prob- 
lem combined with the fact that Pg(x) = (—l) Pe(—x) 
the anion densities follow as p_^(r) = p + j(r)(—l) for 
the Janus dipole. The system is in the sufficiently di- 
lute and weak-coupling regime for the ion-ion interactions 
(a m 0.06), which explains the good agreement between 
both methods. Here we use the association-parameter 
a e [0, 1] introduced in Ref. [H, which gives the equilib- 
rium fraction of the available ions in the electrolyte that 
have formed pairs, to quantify the extent to which strong- 
coupling effects occur. The definition of a in terms of our 
variables reads 



a = 1 



K = 2 



2Kp, 



l - (v/rn^-i); 



dr r exp 



2A, 



(16) 



(17) 



where the fitted value for p s is used. Here K is the 
equilibrium constant for the formation of Bjerrum pairs: 
dipole-like clusters of two oppositely charged ions that 
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Figure 1. A comparison between MC and PB results showing 
the ion densities around a Janus dipole for N± = 300, Ab = d, 
and a — IOAb- In (a) the contour plot shows the net charge 
density p+(r) — P-(r): the MC result (left) and the nonlinear 
PB solution (right). In (b) the same data is represented us- 
ing a Fourier-Legendre mode expansion of the charge density 
p±,e(r), for I = 0, 1, 2, and 3. Note that the negative modes 
in (b) can in this case be mapped onto the positive modes by 
multiplication with ( — 1)*. 

are closely bound due to the strong interaction en- 
erg y 35 i 36 . In Ref. HH it was shown that a < 0.5 im- 
plies that strong-coupling effects are not relevant, i.e., 
the lower the value of a the lower the concentration of 
Bjcrrum pairs (for a = 1 all ions have formed pairs and 
higher order clusters). 

The local dimensionlcss charge density for an equiva- 
lent, homogeneously charged colloid is y = ZAb/(ko 2 ) = 
47tctAb/k ~ 5.2 for the parameters of Fig.[TJ The param- 
eter y can be used to estimate the level of nonlinearity 
in the system and since y ps 5.2 exceeds unity, mode cou- 
pling occur o 45 i 68 . We will come back to this parameter in 
the context of heterogeneous surface charges later. For 
the Janus dipole of Fig. Q] we can clearly observe a non- 
linear effect, namely lim r ^ a p±,o(r) ^ p s , despite the fact 
that the charge on the colloid has no intrinsic monopole 
component. Nonvanishing quadrupole (£ = 2) modes are 
also induced by mode coupling (nonlinearity)^. 

In order to quantify the difference between results ob- 
tained by MC simulations and by PB theory we compare 
the difference in the distribution of ions in the double 
layer directly for each mode. To that end we introduce 
the so-called difference functions fg, which can be ap- 
plied to a general Janus particle with Qtj unit charges on 
the upper hemisphere and Ql unit charges on the lower 
hemisphere, respectively as 

4tt f L/2 
ft = iTTTTlTn / drr \P*Mc{r) - pe,PB(r)\ , 

(18) 

where the ionic charge modes are defined by pt(r) = 
p +t g(r) — p- : i{r) and where we use the labels MC and 
PB to indicate the origin of the respective profiles. Equa- 
tion (jlijp has the property that all fe are when the two 
profiles are exactly the same and that at least one fe>0 



when they are not. Because we compare results for the 
cubic geometry of the simulation box to the spherical ge- 
ometry of the WS cell in PB theory, the upper integration 
boundary is set to L/2 < R. In principle the difference in 
shape and associated boundary conditions imply that we 
compare a simple-cubic crystal of colloids with a liquid of 
colloids at the same volume fraction. However, due to the 
separation of the particles and the level of ionic screening 
the results are virtually independent of the shape of the 
volume when we compare up to r = L/2. This is the 
reason why we only consider systems with added salt. 

In Eq. (JT5J) the functions fe are 'normalized' by |Qul + 
\Ql,\ such that for a homogeneously charged particle 
fo = 2, for the worst-case scenario of full discrepancy 
between the MC and PB results. Because the counter 
charge in the double layer should compensate for the net 
charge on a colloid, each |po( r )| separately contributes at 
most | Qui + |Ql|i which explains the normalization. An 
example of a significant mismatch between the MC and 
PB results is found deep in the strong-coupling regime 
where the MC method predicts a total condensation of 
counter ions in a very small region close to the surface, 
whilst PB theory predicts that the counter charge is lo- 
cated in a diffuse layer around the colloid of significant 
width. For higher order modes the value of fe is bounded, 
but the range is not necessarily [0,2]. To get a feeling 
for the order of magnitude of fe in the case of a good 
agreement between PB and MC results, we mention that 
fo ps 0.003, /i ps 0.073, f 2 ps 0.006, and f 3 ps 0.057 
for the parameter set shown in Fig. [I] whilst Fig. [2ji 
shows profiles with f ps 0.99, /i ps 0.027, f 2 ss 0.034, 
and / 3 ps 0.031. Note that although f in Fig. [3H con- 
firms the huge mismatch between PB and MC results 
for the t — mode, the fe with I > are still small. 
This suggests that the spherical symmetry of the problem 
is only slightly broken within MC simulations and thus 
that strong coupling does not induce significant multi- 
polar charge distributions inside the screening cloud of a 
homogeneously charged particle. 



B. Homogeneously Charged Spherical Particles 

To prove that the difference functions introduced in 
Eq. (fT8|) give a useful description of the deviation between 
the MC and PB results, we investigate the monopole de- 
viation fo for homogeneously charged spherical particles, 
in Fig. [2^,. We compare the behaviour of fo to field- 
theoretical prediction a 57 ' 58 , which are also aimed at es- 
tablishing a range of validity for PB theory, and show 
that there is a good correspondence between the two 
ranges. 

In Refs. [57] andl58l the parameter regimes are investi- 
gated, for which various theoretical approximations give 
trustworthy results for the effective ion distribution of a 
homogeneously charged flat surface. For these flat sur- 
faces parameter space is partitioned into three pieces, see 
Fig. [^Ji. (I) A region where the Dcbyc-Hiickel (DH) ap- 
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Figure 2. A comparison of the data obtained by PB theory 
and by MC simulations of homogeneously charged spheres 
according to the difference function fo of Eq. (|18[) . which 
quantifies the deviation in the distribution of charge in the 
ionic double layer. We show fo as a function of /-c/i, the ra- 
tio of the Gouy-Chapmann and the Debye length, and E, 
the strong-coupling parameter, for several of the systems we 
studied. The field-theoretical prediction of Refs. [H?] and [H 
for homogeneously charged flat surfaces partitions parameter 
space into three regimes, as is indicated by the continuous and 
the dashed line. The Debye- Hiickel (DH), Poisson-Boltzmann 
(PB), and Strong- Coupling (SC) approximations should be 
used to obtain acceptable results in the respective domains. 
The value of the PB-SC divide is denoted by H* « 10. (b-d) 
Three samples of the ion profiles that are obtained by PB the- 
ory (full curves) and MC simulations (dashed curves), showing 
cation and anion densities. Here fo = 0.03, 0.39, and 0.99. 
The discontinuity in the first derivative of the MC p+,o profile 
in (d) is caused by the positive ions condensing on the surface 
of the colloid in combination with the binning procedure we 
applied to average the ion densities. A similar discontinuity 
is present in the p_ : o results, although this is not visible on 
the scale of this plot. 



proximation^ 9 - can be used, the screening is linear. (II) A 
region where the charge of the surfaces becomes higher, 
the nonlinear PB equatio n 18 ! 19 has to be solved in order 
to determine the effective electrostatic interactions. (Ill) 
A region in which the ion-ion correlations close to the sur- 
face require the use of Strong-Coupling (SC) thcory22r— . 



In Fig.[2]thc parameters k/i = 2/y and S = (y/2)kAb rep- 
resent parameter space in a 'field-theoretical language' 
(see Ref. [El), with y 'our' local dimensionless charge den- 
sity [as defined just below Eq. (fTT))] and \x the Gouy- 
Chapmann length. The Gouy-Chapmann length must 
not to be confused with the ionic chemical potential, 
which is also usually denoted by [i. S is a measure for the 
particle's surface charge, expressed in terms of the Bjcr- 
rum length. PB theory produces satisfactory results for 
5 < S*, with S* a cut-off value: Ref. HH sets the value 
of S* 10 for the transition between the PB and SC 
regimes. The value of Kfj, is of minor importance when 
H < S* as PB theory can be straightforwardly applied to 
the DH region in the low-charge limit. 

By comparing the MC and PB ion profiles for the 99 
systems containing a homogeneously charged sphere that 
we have investigated, we found /g = 0.1 to be a good 
boundary value for the regime (/g < /g) in which PB- 
theory accurately describes the ion density profiles. Fig- 
ures EJ^b-d) show three example ion profiles to illustrate 
the possible level of deviation corresponding to a particu- 
lar fo value: /g = 0.05, 0.39, and 0.99, respectively. Note 
that our choice of /g = 0.1 is meant to imply that below 
this value PB theory gives a good description, however, 
for /g > /g PB theory may give a reasonably accurate 
description (see Fig. [2}:), but this is not a given and SC 
theory may be required to give a good description. Using 
our /g values as a function of y and S, see Fig. [2k, we 
locate the PB-SC divide at S* ss 1. Minor changes in 
the value of /g do not significantly change the location 
of the PB-SC transition in parameter space. However, 
since what is considered an unacceptable level of the dis- 
crepancy between PB and MC results is dependent on 
the quantities/behaviour we are interested in, there is a 
degree of arbitrariness to our result. Nevertheless, our 
approach to this problem and our choice for /g appears 
justified since we obtain a similar partitioning of parame- 
ter space as was found in Ref.[58|. This was to be expected 
for a homogeneously charged sphere, since there is only 
a geometrical difference with respect to a homogeneously 
charged plate, which for sufficiently large spheres can be 
considered small close to the sphere's surface. Our results 
show that even for relatively small spheres (compared to 
the size of the ions) there is qualitative agreement. 

For completeness we comment on the accuracy of our 
MC result deep inside the strong-coupling regime. The 
MC results show that a layer of oppositely charged ions 
can form on the surface of the charged particle. The 
interaction between the charges (sites and ions) is such 
that the free ions effectively condense on the particle, 
see Refs. HH, [70l - t73l for a more comprehensive account 
of this phenomenology The ions in the electrolyte ex- 
perience similarly strong interactions and form Bjerrum 
pairs. Since we only consider single particle MC trial 
moves, the formation of Bjerrum pairs interferes with the 
exploration of phase space in the strong-coupling limit. 
The clusters hardly move, because most single particle 
moves that would break up a cluster are rejected based on 
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the energy difference. This results in an ill-converged en- 
semble average, when the Bjerrum-pair concentration is 
high (a > 0.5) 3 - 6 -. The problem can be overcome by intro- 
ducing cluster and association-dissociation moves for the 
Bjerrum pairs to obtain a more efficient samplin g 35 ' 36 . 
However, we do not believe that ion condensation and 
Bjerrum-pair formation will influence our result with re- 
gard to the location of 5*, since these effects only start 
to play a role for 3 ^> S*, since then a > 0.5. 



C. Janus-Dipole Charge Distributions 

1. Comparison of the Monte Carlo and 
Poisson-Boltzmann Results 

In order to investigate the range of applicability of the 
(nonlinear) Poisson-Boltzmann approximation, we apply 
our method of comparison from Section ()III A[) to higher 
order Fouricr-Legendre (FL) modes of the ion density for 
the case of a Janus dipole. Figure [3J shows the deviation 
parameter ft for I = 1,3 and the large set of param- 
eters we studied. Note that the Janus dipole bears no 
even multipole moments \p-.i(r) = p + t i{r){— l) e ] due to 
its antisymmetric charge distribution. To app ly a rep- 
resentation similar to the one used in Ref. [57| for Janus 
particles, we introduce the following modified dimension- 
less parameters: y s = 2/(k/x s ) = (|Qu| + IQlDAbAko 2 ) 
and = (ys/2)kAb. The sum of the absolute value 
of the charge on each hemisphere is used, rather than 
the total charge (which would be zero in the case of a 
Janus dipole). For pure monopoles ys and 2^ reduce to 
the original parameters y and 2. We prefer to express 
our results in terms of the dimensionless (absolute) local 
charge density ys rather than in terms of k/j, = 2/ys, 
since the former is a more natural quantity for PB the- 
ory of colloid systems. The use of k/j, in Fig. [2] was to 
illustrate the correspondence between our results and the 
partitioning given in Ref. [58l 

For the dipole term (£ = 1) we observe trends in the 
value of fi (Fig. [3]) as a function of y and 2 similar to 
those observed for the value of f of a homogeneously 
charged sphere (Fig. [5^). The onset of a strong difference 
in the correspondence between the two results for the 
dipole mode occurs at /* ps 0.1. For the octupole (£ = 3) 
term, the crossover value /jT for an appreciable level of 
deviation appears to be slightly larger than 0.1, but on 
the strength of our results it is difficult to state this with 
certainty. 

Based on Fig. [31 a regime can be distinguished for the 
leading dipole term where PB theory yields accurate re- 
sults for the charge profiles in the electric double layer 
(He < Sg ps 1). For the I = 1 through £ = 5 modes 
(£ = 5 not shown here) the correspondence between MC 
and PB results is also sufficient when 2^ < S E . The 
modified parameter 2^ therefore appears useful to de- 
scribe parameter space for dipolar Janus particles with 
regards to quantifying the region where PB theory can 
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Figure 3. The deviation fe in the double layer, determined us- 
ing the MC simulations and PB theory, for a Janus dipole as 
a function of the modified charge density ya and the modified 
strong-coupling parameter He. The subgraphs show the re- 
sults for the first two odd FL modes (I — 1, 3), corresponding 
to the dipole and octupole term, respectively. 



be used to describe the system. 



2. Multipole-Expanded DLVO Approximation for the 
Janus Dipole 



Equations describing the electrostatic pair-interaction 
of (spherical) colloidal particles with a substantial dipolar 
contribution to the surface charge are derived in the Ap- 
pendix. This derivation is performed within the Poisson- 
Boltzmann approximation, and the resulting equations 
can therefore be considered as an extension of DLVO 
theory towards inhomogeneously charged particles. The 
monopolc-monopolc, monopolc-dipolc, and dipolc-dipolc 
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interactions are given respectively by 

,MM /r> ^ _ i eXp( — fjRjj) ^ Y ryY 



/3V^ M (R ij ) = A E 
exp(-Kftj) 

B R2 

/3V° D (Rjj) = 
exp(— nRi 



Ri 



((p7 ■ R tJ )Zj - (pj ■ Ay)Z? 



(19a) 



(19b) 



i? 3 . 



(' 



l + /c.Ry)(p7 -pj) 



(3 + + («^) 2 )(pf ■ iSy)(pJ ■ Rij) 



(19c) 



with Ry the distance vector between particle i and j, 



| |, and Rij — R^- /Rij. We also introduce 



the ' Yukawa- monopoles' and 'Yukawa-dipoles' p^, 
which, for Janus spheres of radius a, are given by 



Pi 



(Qu 
(Qu 



Qi 



Qi 



cxp(«;a) 



(20a) 



3a exp(Ka)Ai 



4(2 + e c /e)(l + Ka) + (Ka) 2 ' 



(20b) 



with iii is the particle's symmetry axis, which points to 
the northern hemisphere, Qu the total charge on the up- 
per hemisphere, Ql the total charge on the lower hemi- 
sphere, and e c /e the ratio between the relative dielectric 
constant of the colloidal particle and that of the sur- 
rounding medium, which we choose 1 unless stated dif- 
ferently. However, Eqs. (fT9|) and ([20)) arc, as is typical 
for DLVO theory, only valid for sufficiently small charges, 
since linearised PB theory (Debye-Hiickel approxima- 
tion) was employed in the derivation of these equations. 
Fortunately, charge renormalisation has proven to be a 
useful tool in broadening the range of applicability to- 
wards particles with a higher charg e 19 ' 24 . This renor- 
malisation procedure was extended towards dipoles and 
higher multipoles in Ref. l45l 



hemispherical charge distribution than for the other two 
distributions we considered. This is because in the MC 
simulations we used 50 divalent charge sites on the upper 
hemisphere instead of monovalent sites. Our results for 
the hemisphere are therefore less convincing than for the 
other two charge distributions, but our preliminary indi- 
cation is that the range in which the PB approximation 
is valid is roughly the same as that for the monopole and 
Janus dipole. 



2. Multipole-Expanded DLVO Approximation for the 
Hemispherical Charge Distribution 

As mentioned in the previous section, one cannot 
always resort to standard DLVO theory for particles 
with an anisotropic surface charge, because dipole, 
quadrupole, or higher order multipole charge distribu- 
tions are relevant. Consequently, two-body interactions 
of the form monopole-dipole, dipole-dipole, monopole- 
quadrupole, etc., should be considered as well. 





D. Hemispherical Charge Distributions 

1. Comparison of the Monte Carlo and 
Poisson-Boltzmann Results 

For the hemispherical charge distribution, we also com- 
pared PB results with those of MC simulations, in order 
to determine the regime of validity of the PB-multipolc 
expansion. We again considered 99 systems and found 
that for the even modes the level of deviation ff m 0.1 
sets a rough upper bound to the applicability of PB the- 
ory. For the odd modes the correspondence between the 
MC and PB results seems to hold for slightly higher val- 
ues of the deviation parameter: f£ ps 0.25. Using these 
two values of /* we can roughly locate the range of va- 
lidity of the PB result in the region Be < 1. The effects 
of strong coupling are however far more apparent for our 



Figure 4. (a) A sketch of two interacting particles with radius 
a, both having charge on only one hemisphere. To eliminate 
the dipole moment we relocate the centre of the charge dis- 
tribution at a distance b (along the particle's rotational sym- 
metry axis) from the geometrical (hard-core) centre of the 
particle. This choice results in a distance Rij between the 
charge distributions. Graphs (b) and (c) show the suggested 
value of b and the Yukawa weight factor C(kcl), respectively, 
as a function of the colloid radius a in terms of the Debye 
screening length for a wide range in e c /e, the ratio be- 
tween the relative dielectric constant of the particle and that 
of the surrounding medium. In (c), we indicate C = 1, which 
is the weight in case of the regular DLVO equation, using a 
dashed line. 

Figure 0^ shows two spherical particles with a positive 
hemispherical charge distribution, i.e., one side is charged 
the other is uncharged. Instead of choosing the centre of 
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the charge distribution to coincide with the geometrical 
centre of a spherical particle one is free to place this point 
anywhere inside the particle. The most natural location 
is the point for which the (Yukawa- )dipole moment van- 
ishes and thereby all dipole interactions. With this choice 
the electrostatic two-body interaction in terms of only a 
monopolc-monopolc term, is expected to be maximally 
accurate. This point is located on the rotational sym- 
metry axis of the hemisphcrically charged colloid. For 
sufficiently large interparticle distances we obtain the fol- 
lowing interaction potential for the shifted-monopole ap- 
proximation: 

-^(%) * (z-^-C( K a)) 2 AB ex P(-*%) ; (21) 

with the shifted centre-to-centre distance, Z the par- 
ticle charge, and C{na) a renormalisation factor. This 
renormalisation factor depends on the ratio of the parti- 
cle radius a and the Debye screening length only and 
C = 1 in case of regular DLVO theory. We can numeri- 
cally determine the distance b from the particle's geomet- 
rical centre, where we should place the shifted monopolc 
such that the dipole term vanishes. Figure Hp shows the 
ratio b/a as a function na for several values of e c /e. Note 
that b/a goes to 1/2 for small na, which is the (Coulomb) 
result that is obtained for unscreened particles. We 
show the corresponding renormalisation factors C{na) in 
Fig. EI C )- We find C(na) < 1, because the charge is 
located closer to the (new) centre of the charge distri- 
bution than for homogeneously charged particles. Note 
that the interaction potential is independent of the orien- 
tation of the particles, when rotated around the shifted 
(monopolc-charge) centre. The monopolc-monopole ap- 
proximation therefore does not capture the full interac- 
tion between the two hemispherically charged colloids. 
To capture the orientational dependence higher order 
modes (octopole and higher) are required. However, for a 
monopole-monopole only approximation Eq. (|2ip is max- 
imally accurate by design. In the Appendix we explain 
how to calculate and, thereby, also how to minimize the 
Yukawa dipole moment for any type of charge distribu- 
tion and hard-core shape. 



IV. CONCLUSION AND OUTLOOK 

In this paper we studied the range in parameter space 
for which the nonlinear Poisson-Boltzmann (PB) theory 
accurately describes the behaviour of the ions around a 
Janus charge-patterned spherical colloid in a 1:1 elec- 
trolyte. We used primitive-model Monte Carlo (MC) 
simulations to establish the ion density around such a 
charged particle for a huge set of parameters. By also 
computing the ion density for the same parameters and 
comparing the two results, we were able to establish a 
regime in which this PB theory gives a good approxima- 
tion for the ion distribution. This comparison is based 



on Fourier-Legendre (FL) decomposition of the MC ion 
density to determine the contribution of the monopolc, 
dipole, quadrupole, . . . charge terms. The theoretical ap- 
proach also relies on FL decomposition and this enables 
us to quantify the differences on a mode-by-mode basis. 

For a homogeneously charged sphere we compared our 
range of validity for PB theory to the range found in 
Refs. H3 and [Hf for a system of homogeneously charged 
flat plates in an electrolyte. There is a remarkable corre- 
spondence between the two ranges, especially considering 
the small size of the colloids that we studied in relation 
to the size of the ions. For such small spheres a greater 
deviation with respect to the results of a flat-plate calcu- 
lation could reasonably be expected. We were also able 
to show that the range in which the PB results accu- 
rately describe the ion density around a spherical Janus- 
dipole is similar to that found for the homogeneously 
charged sphere. For particles with only one (homoge- 
neously) charged hemisphere, there is an indication that 
the regime in which PB theory can be applied matches 
the regime found for the two other particles. 

In the PB-regime that we obtained, we can use sim- 
ple (multipole-expanded) DLVO-likc equations, which we 
derived in this paper, to describe the interactions be- 
tween two particles with a Janus-type charge. We gave 
explicit expressions for the monopole and dipole inter- 
actions, since these terms are typically dominant for 
Janus particles. These electrostatic interactions resem- 
ble well-known Yukawa interactions, and reduce to these 
in the homogeneous charge limit. For the Janus-dipole 
obtaining the (multipole-expanded) DLVO expression is 
relatively simple. To accurately model a hemispherical 
charge distribution using only a monopole-term is a lit- 
tle more complicated. The key step proved to be shifting 
the centre of the charge distribution from the centre of 
the particle towards the charged hemisphere in order to 
eliminate the Yukawa-dipole contribution. 

Our analysis forms a basis of a good understanding 
of the range in parameter space for which the PB ap- 
proximation can be applied to describe the behaviour of 
heterogeneously charged colloids. This is, for instance, 
relevant to the study of such particles using simulations, 
where PB-theory-based effective interactions can be used 
to study the phase behaviour of such particles in the 
right regime. Note that we only considered equilibrium 
ion density profiles of stationary colloids. The rotational 
movement of mobile charge-patterned colloids can occur 
on time scales that would lead to an out-of-equilibrium 
double layer. What effect the out-of-equilibrium ion den- 
sity would have on the screening of the particle and how 
such effects should be incorporated into effective inter- 
action potentials used in simulations, is left for future 
investigation. 
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APPENDIX: JANUS PARTICLES AND THE 
DEBYE-HUCKEL APPROXIMATION - DLVO THEORY 
FOR PATCHY COLLOIDS 

In this appendix we derive the (multipole-expandcd) 
DLVO-like equations for general charge distributions and 
we show how this leads to Eqs. (fTi?)) - ([2"0")l . Before that we 
set the stage by first examining the quality of the dipole- 
only approximation for a Janus-type charge distribution. 




Figure 5. A contour plot of the net charge density p+(r) — 
p- (r) around an antisymmetric Janus particle for the param- 
eters Z = 10, N± = 425, Ab = d, and a — IOAb, showing the 
profile that follows from a mode expansion up to I — 6 on the 
left, on the right only the dipole mode is plotted. 

We have shown that the mode expansion up to £ = 6 
gives good agreement between MC simulations and PB 
theory regarding the ion profiles in a well-defined regime 
of parameters. For Janus particles in general the most 
dominant multipolcs are the monopole and/or the dipole; 
note, however, that for 'pure' Janus dipolcs the monopole 
term vanishes. To show this Fig. [5] plots the PB result 
for the local ionic charge densities around such a purely 
dipolar Janus particle, using a smaller charge (Z = 10) 
than elsewhere in this paper. The number of ions added 
to the system, N± = 425, corresponds to na w 2.8. We 
find that ys = 0.36 and 5s = 0.05 for this system, which 
implies that we are within the regime where PB-theory is 
applicable according to our analysis. We are aware that 
a multipole expansion of Yukawa-like interaction poten- 
tials does not necessarily converge in genera l 74 ' 75 . In this 
particular case, however, we catch most of the physics of 
the interacting Janus particles, by treating the particle 
as a 'pure' dipole (without higher order modes). To il- 
lustrate this, the left side of Fig. [5] shows the ion charge 
density up to I = 6, which we proved to be sufficient 



for good correspondence with MC results in the regime 
where PB-theory is applicable. The right side shows the 
dipole {1=1) mode only. The differences are therefore 
due to the missing I = 3 and t = 5 terms. The dipole 
approximation overestimates the electrostatic potential 
in the axial direction, whilst it underestimates the elec- 
trostatic potential in the perpendicular direction. 

With this in mind we explain the way to extend the 
applicability of DLVO theory to anisotropically charged 
particles. The result of setting up this theory allows us 
us to find explicit equations for the monopole-dipole and 
dipole-dipole interaction potential between Janus parti- 
cles as a function of their orientation. These expres- 
sion are similar to the well-known expressions for the 
interaction of unscreened dipolcs. We begin by consid- 
ering the effective electrostatic energy of an extended 
charge configuration eq(r) in a 1:1 electrolyte with bulk 
concentration 2p s - we do not incorporate hard-core ef- 
fects at present. By using the electrostatic energy from 
Coulomb's law combined with the ideal-gas entropy for 
the monovalent ions, we find that the grand potential of 
the charge configuration inside a two-component mono- 
valent ion mixture inside a solvent is given by 

+ i J dr (p+(r) - p_( T ) + <z(r))0(r), (22) 

where p±(r) are the ion densities and where the dimen- 
sionless electrostatic potential is 

W-^/d ^W-^W + rtO , (23) 

Note that we applied the Debye-Hiickel approximation; 
the 'usual' entropic term of the ions is linearised via 
P±(r)(log(p±(r)/p s ) - 1) ~ P±(v)(p±(r)/ps ~ 2), imply- 
ing that we assume the ion densities not to vary much 
from the bulk value p s . For a more detailed derivation of 
Eq. ([22]) see for example Ref. [(ill 

We consider a system that consists of a collection of 
M (fixed) charges located at r c j, with i = 1, . . . , M an 
index. These charges should not be considered as point 
charges, but rather as localized charge distributions <?i(r) 
inside associated volumes Vi, which are non overlapping 
and centred at r c j, as is sketched in Fig. El Each qi(r) is 
only nonzero inside Vi, and the total (non- ionic) charge 
density is hence written as q(r) = £ i=1 9i( r )- 

Since wc have not included any hard-cores yet, the ion 
densities follow immediately from setting the functional 
derivative of Eq. (f2"2"j) w.r.t. the ion densities to zero, 
SH/6p ± (r) = 0, yielding p±(r) = p s {l T &('))> 
with 

*(r)=A B / dr-^rO ^T^r^ ' ™ 
Jv, \ r - r l 

where k 2 = 8it\bp s . By assuming that the charges qi(v) 
within the volumes Vi have a fixed position, Eq. (|22[) can 
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Figure 6. A sketch of two interacting charge distributions 
(grey) 5i(r) and ?j(r), inside the enclosing volumes Vi and 
Vj , respectively. These volumes are rod- like in this particular 
example to resemble rods with one charged 'head' and have 
boundaries that are indicated by dashed instead of solid lines 
since we do not consider hard cores while calculating the ion 
densities connected to these charge distributions. The centres 
of the charge distributions can be chosen arbitrarily and are 
indicated by and Yj. 



Eq (I22j) may therefore be rewritten as 

PH=f dr 5>°^(r)(^M-2 
Jv out ^± V Ps 



dr -( P T(r) - P°- Ut (r) + g(r))0(r), (26) 



with 

4>(r) = A B J dr — — , (27) 

and the 'new' colloidal charge distributions 
q(r)=q(r) + p i ?(r)-p i «(r) 

M 

= E*W- (28) 



be written, up to a self-energy constant, as 



pH = A B V 

i<j JJ VuV, 



A A ' ( \ I >\ eX P(~ K \ r ' - r \) 

drdr ft(r)^(r) ^— ^ . 

(25) 

Equations (f2~l| and (|25p were derived using the Taylor- 
expanded (quadratic) Hamiltonian ([22]) and therefore 
give appropriate results only in case the dimension- 
less electrostatic potential </>(r) remains sufficiently small 
w.r.t. unity. 

We will now use the result of Eq. (f2"5j) to obtain an 
analogous theory for charge distributions with associated 
hard-core volumes. This is done by 'freezing' the ionic 
charge profiles inside the volumes Vi and adding these to 
the charge configurations <?i(r), such that the new charge 
distributions (r) are found. Effectively we compensate 
for the fact that ions in the Yukawa approximation can 
penetrate the hard particle. We thus consider a new com- 
bined charge distribution ft(r) that is exactly the distri- 
bution we are interested in, by compensating the fixed 
charge for the ion profiles it induces. 

Starting with the obtained ion densities for the system 
without hard cores, we split the entire system volume V 
into Vout, consisting of all points r outside the volumes 
Vi for all i, and the volume V; n of points that are inside 
one of these volumes. Note that V n is the complement 
of V ou t- The ion densities are also split into /9j_ ut (r) and 
pf(r) such that /4 ut (r) + p^(r) = p±(r) and p| ut (r) = 
for all r inside V m: and p±(r) = for all r inside V ut- 



for which g»(r) = qi(r) + p+(r) — /c™(r) when r inside Vi. 

Within the DLVO-approximation, the net ionic charge 
density p™(r) — p™(r) in volume Vi is only induced by 
the fixed charges <?i(r) in that volume itself (large inter- 
particle distances). One therefore finds for r in Vi that 

<ft(r) « g-i(r) - 2p s 0.;(r) 

= gi (r)-2ABp s /dr- gi (rO eXp( r; |r/ r r|) - 



v, 



(29) 



The second line in Eq. (|26| becomes a constant that may 
regarded as a self-energy term and therefore can be ig- 
nored to yield 



[3H= I dr V>±(r) 

./Vout -L 



a=± 



P±( r ) 



1 



+ - J dr (p + (r) - P -(v) + ?(r))0(r), (30) 

with the restriction that p±(r) = if r is inside V m . Com- 
paring Eq. (|22jl with Eq. (|30|) . the effective interaction 
Hamiltonian between hard particles with charge densi- 
ties <?i(r) can be recognized in the latter. The effective 
interaction energy for such a system can thus be obtained 
from Eq. with gi (r) the solution of Eq. (|2"5]). with 
5i(r) the actual charge density of interest. For particles 
with hard-core volume Vi and corresponding charge den- 
sities q~i(r) it is necessary to first solve qi(r) from Eq. (|29|) . 
By finding qi(r) we obtain the charge distribution that 
compensates for the self- induced ion densities, since these 
are required in to determine the effective interaction in 
Eq. (Ell). 

At this point we have found a procedure to calculate 
the electrostatic interactions between inhomogencously 
charged colloidal particles with hard cores, within the 
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Vi 



Figure 7. A sketch of two interacting charge distributions qt 
and qj with hard-core volumes Vi and Vj, respectively, sepa- 
rated by a distance Rij . We also show the vectors s; and Sj , 
which have their origins at r c i and r CJ -, respectively. 



DLVO approximation that we discussed earlier. In order 
to isolate specific multipole interactions between colloids, 
Eq. (f2"5")) is rewritten as 

PH = A B y~] II dsjdsj \qi(r c i + Si)qj(r cj +Sj) 

i<j JJ Vi,V,3 I 

exp(-/«|Ry + - Sj|) 



R, 



(31) 



with = r C j — r ci the colloid-colloid distance and 
Sj = r — r ci the coordinate relative to r c?; . Both arc shown 
in Fig. [7] The Yukawa interaction can then be expanded 
into spherical harmonics around these centres^, which 
are r C i and r C j. For this we must assume that the dis- 
tance from any point in Vi to its centre r C i is less than the 
distance to any other centre r C j . This however automati- 
cally holds for equi-sized spherical volumes if one chooses 
the centres r C i exactly in the middle of the spheres. 

Here, the monopole and dipole terms are of main inter- 
est. Eq. (|31[) can therefore be written, up to a monopolc- 
monopolc, monopole-dipolc, and dipolc-dipolc interac- 
tion term as 

// X V.^ M (%) + V^";R„! + V£ D (Ry)- (32) 

i<j 

The resulting interaction terms are given by Eq. (|19l) , in 
which the 'Yukawa monopole' and 'Yukawa dipole' are 
now (for the general charge distribution) 



Zj = / ds,; qi(r ct + s,) 
JVi 



sinh KSi 



h S , 



(33a) 




(33b) 



with Si = |s,|. Note that Eqs. (|33a[) and (|33b[) reduce to 
the well-known expressions for Coulomb monopole and 



dipole in the limit na 1 0, sec Ref. [77|. Also the in- 
teraction terms (fTi?|) give the well-known result for pure 
Coulomb systems in this limit. Finally, the electrostatic 
potential around charge distribution i can be expanded 
as 

(pi{r ci + Si) = Z { A B 



. Y - n > exp( — KSi) N 
+ (pj ■ S, A B — ^ "(I + KSi) 

+ 0{£>2), 



(34) 



in which quadrupole and higher order multipoles are not 
included. 

As an example we show that the interaction between 
particles with a homogeneous charge distribution is in 
agreement with DLVO theory. By considering a spher- 
ical colloid with a hard core radius a and a charge Z 
distributed homogeneously over its surface and choosing 
r c i = for convenience, the charge distribution is given 
by q~i — Z/(4-Ka 2 )S(si — a). It can be shown that Eq. ([29]) 
is solved by 



9i(si) 



■ z 





S( Si - a )+ ^l 2 ?* 

V 1 1 ' a{l+K,a) 



if Sj 
if s. 



<a, 
> a. 



(35) 



Note that the 'added' charge density inside the colloids 
has the same sign as the surface charge and exactly can- 
cels the ionic charge due screening. From Eq. (|33a[) one 
obtains Zj = Zexp(Ka)/(l + na) and this gives the 
DLVO result as Eq. (|19aj) shows. 

The Yukawa monopole and the Yukawa dipole can also 
be extracted from the solution for the electrostatic poten- 
tial outside a spherical particle, using Eq. (|34[) . without 
calculating the charge density qi(r) explicitly. Namely, 
we can often solve the (linearized) Poisson-Boltzmann 
equation around a single particle mode-by-mode, and 
then read the multipole moments from the mode am- 
plitudes in the final solution. As an example we con- 
sider spherical particles with a fixed surface charge dis- 
tribution that is rotationally symmetric around the unit 
vector hi, e.g., a Janus particle. The electrostatic po- 
tential in the vicinity of the single particle can be ex- 
panded as 4>(s u Xi) — XX (f>i(si)Pi(xi), withal = s r nj, 
and Pe(xi) the £'th order Legendre polynomial. Now 
( = and 1=1 correspond to the monopole and the 
dipole contribution to the electrostatic potential, respec- 
tively. The multipole mode functions, which solve the 
linearised PB equation, behave as (j>e(si) ~ s\ for r < a 
and 4>e(si) ~ ki(nSi) for r > a, with ki the i'th modi- 
fied spherical Bcsscl function. By applying the boundary 
condition 

lim <p'(si,Xi) = lim <p'(si,Xi) - 47rA B cr(xi), (36) 

Si-la Si^a 

with the prime (') denoting the radial derivative (w.r.t. 
Si), and 17(2;,) the surface-charge density, a solution to the 
electrostatic potential in terms of a multipole expansion 
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is obtained. Since the modified spherical Bcssel functions 
can also be recognized in Eq. (jM]) . one is able to derive 
that the monopole and the dipole moment should relate 
to the expansion of the surface-charge density, a(xi) = 



7 y A 2 exp(Ka) 
Zi = 47ra ctq 



pj = 47ra 3 CTi 



1 + Kfl 

exp(rea) 
3 + 3na + (na) : 



■n, 



(37) 
(38) 



In the case of Janus particles with charge densities <j\j = 
Qu/(47ra 2 ) and ctl = Qh/(^a 2 ) on the upper and lower 
hemisphere respectively, one finds Oq = ^(cy + ol) and 

Throughout this paper we assumed that the dielec- 
tric constant of the particles matched the dielectric con- 
stant of the solvent. However, Eq. (|38|) can easily be ex- 
tended to particles that have a dielectric mismatch with 
the solvent. We modify Eq. (j38]l by including a dielectric 
jump in the boundary condition [Eq. (|36[) ]. such that the 
Yukawa multipoles for these particles can be obtained as 
well. As is known from the DLVO equation, we find that 
the Yukawa monopole is unaffected by the value of the 
dielectric constant inside the particle, whilst the Yukawa 
dipole changes. Equation (|38|) becomes 



= 47ra 3 <7i 



cxp(Ka) 



{2 + e c /e)(l + na) + (na) 2 



(39) 



with e c /e the ratio of the relative dielectric constant in- 
side and outside the particles. Note that the dipole 
moment becomes very small if the interior of the col- 
loid is very polar w.r.t. the surrounding medium, e.g., 
(Pickering) emulsions of water in oil. Consequently the 
monopolc-dipole and dipolc-dipole interactions will be 
negligible for these systems, even in case of an asymmet- 
ric charge distribution. However, if the dielectric con- 
stant of the colloid is comparable or smaller than the 
dielectric constant of the solvent a significant dipolar in- 
teraction may arise. 
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